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Drude weight {D) is a useful measure to distinguish a metal from an insulator. However, D has not been justifiably 
estimated by the variation theory for long, since Millis and Coppersmith [Phys. Rev. B 43 (1991) 13770] pointed out that 
a variational wave function 'Pq, which includes the key ingredient (doublon-holon binding effect) for a Mott transition, 
yields a positive D (namely metallic) even in the Mott-insulating regime. We argue that, to obtain a correct D, an 
imaginary part must exist in the wave function. By introducing a configuration-dependent phase factor Vq to 'Pg, Mott 
transitions are successfully represented by D (D = 0 for t/ > Uf) for a normal and J-wave pairing states; thereby, the 
problem of Millis and Coppersmith is settled. Generally, Vq plays a pivotal role in describing current-carrying states in 
regimes of Mott physics. On the other hand, we show using a perturbation theory, the one-body (mean-field) part of the 
wave function should be complex for band insulators such as antiferromagnetic states in hypercubic lattices. 
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1. introduction 

The Drude weight D, the coefficient of the DC (oj = 0) 
conductivity, has been considered to be an important measure 
to distinguish a metal (D > 0) from an insulator (D = 0), in 
particular, since Kohn showed that D can be calculated only 
with quantities with respect to the ground state. ^ Actually, D 
is obtained through D = d^E(A)/dA^\A^o (E: total energy) by 
introducing a virtual flux (Peierls phase) A or equivalently by 
twisting the boundary condition. This formalism is not only 
convenient for a variety of approaches^’ ^ but crucial in partic¬ 
ular for the variation theory. 

More than two decades ago, Millis and Coppersmith'^ first 
applied this formalism to variational wave functions for the 
one-dimensional Hubbard model at half filling.^ They used 
wave functions that include binding factors between a dou¬ 
bly occupied site (doublon, D) and an empty site (holon, 
H),^’^ Pq, in addition to the usual onsite correlation (Gutz- 
willer) factor,^ Pq, and showed that this type of wave func¬ 
tions ('Pn = PqPG^N with On being a Fermi sea) are metal¬ 
lic even for a sufficiently large Ujt to be insulating, in the 
sense that D > 0. At that time, it had not been clarified 
yet that a D-H factor Pq brings about a Mott transition, al¬ 
though the Gutzwiller wave function ('Pq = ^g^n) was 
known to be always metallic.^’Consequently, their result 
has caused confusion and misunderstandings—'P n cannot de¬ 
scribe a Mott transition— to subsequent studies using 'Pn 
( and 'Pj = PqPo^d with O^ being a J-wave BCS state). 
Later, it was confirmed^ that 'Pn and similar wave 
functions that have D-H binding effects induce Mott tran¬ 
sitions at (7c ~ W {W\ band width); the behavior in vari¬ 
ous quantities changes with anomalies at (7c, such as dou¬ 
blon density d —an order parameter of Mott transitions—, 
charge-density structure factor A^(q), and momentum distri¬ 
bution function ^(k). Thus, the discrepancy between the be¬ 
haviors of these quantities and D has remained an enigma for 
long years. 

The main purpose of this study is to develop a method for 
calculating the Drude weight appropriately in the variation 
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theory, and settle the above problem. We first argue that a fi¬ 
nite imaginary part is indispensable in the wave function for 
correctly estimating D. 

In the regime of Mott physics ((7 > (7c), we noticed that 
some configuration-dependent phase factor has to be intro¬ 
duced to the ordinary (real) wave functions like 'Pn and 'Pj. 
Analyzing the phase added in hopping between a doublon and 
a holon in the field A, we construct a phase factor for the 
present case Pq, where ^ is a phase parameter to be optimized. 
For Mott insulators, the optimized Pq cancels out the Peierls 
phase by satisfying the relation 6 = A, and the increment of 
energy owing to A is reduced to zero. Thus, D vanishes for 
U > Uc Sit half filling. For (7 < Uc or S > 0 (S: doping rate), 
the Peierls phase is cancelled out only partially, and E{A) re¬ 
mains larger than £"(0); D becomes finite. In calculations for 
A = 0, ^ is optimized at zero {Pq = 1), so that the previous 
results for J, A^(q) and ^(k)^^“^^ remain unchanging for the 
new wave function 'P = Pe^^^ Thus, the long-standing issue 
of Millis-Coppersmith was resolved. 

In fact, we found that this type of configuration-dependent 
phase factors seems generally essential to treat current- 
carrying states appropriately in the regime of Mott phys¬ 
ics ((7 > (7c). It was shown that a similar configuration- 
dependent phase factor is indispensable to correctly represent 
staggered flux states in Hubbard-type models. 

On the other hand, in treating a band insulator for (7 (7c, 

we And Pq is ineffective. For instance, although an antiferro¬ 
magnetic (AF) state for the square-lattice Hubbard model at 
half Ailing is insulating for any Ujt (> 0), the Drude weight 
obtained using 'Paf = PqPqPg^kv becomes finite for small 
values of UIt. It indicates another element is needed for 'Paf- 
We argue using a perturbation theory that the imaginary part 
in the one-body (mean-field) wave function Oaf is vital for 
reducing D in this case. Here, the imaginary part is given by 
the first-order perturbation with respect to A. Thus, the mech¬ 
anism to suppress D is different between Mott and band insu¬ 
lators. 

We also address the effectiveness of Pq for superconducting 
(SC) weight Ds in the attractive Hubbard model. 

This paper is organized as follows: In Sec. 2, we introduce 
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trial wave functions for estimating the Drude weight through 
Kohn’s formula, and mention why 'Fg fails in yielding an ap¬ 
propriate D. In Sec. 3, we give the results of VMC calcula¬ 
tions using for the normal and J-wave pairing states, 

and discuss an improvement of the phase factor. In Sec. 4, we 
consider an AF state as a typical case of band insulators of 
weakly correlations. In Sec. 5, the SC weight is reconsidered 
for the attractive Hubbard model. In Sec. 6, we recapitulate 
the main results in this paper. Preliminary results of this study 
were presented in a proceedings.^^ 


ational parameters are efficiently optimized using the quasi- 
Newton algorithm.^^ In some cases, we adopt the stochastic 
reconfiguration method^^ to reduce statistical fiuctuation. We 
use L X L systems (L = 12-16) with the periodic-antiperiodic 
boundary conditions. Although we implement calculations 
using the square lattice, the properties obtained for the nor¬ 
mal state are qualitatively independent of the form and di¬ 
mensionality of the lattice. In most cases, samples as many 
as 2.4 X 10^ are used to reduce numerical errors, which are 
typically ~ for the total energy. 


2. formulation 

In Sec. 2.1, we introduce the model and method used in this 
study. In Sec. 2.2, we explain ordinary trial wave functions for 
f/ > 0 without a Peierls phase. Trial states for f/ <0 will be 
discussed in Sec. 5. In Sec. 2.3, we show the necessity of a 
complex wave function for a finite Peierls phase, and how to 
construct a trial wave function in regimes of Mott physics. In 
Sec. 2.4, we discuss the relation between Drude (D) and SC 
(Ds) weights in the light of variation theory. 

2.7 Model and method 

To study Mott transitions through the Drude weight, we 
adopt the Hubbard model on the square lattice: 

"H = ^kin + ^int 

= + H.c.) + 1/ y] nj^nji, (1) 

j 

where (/, j) denotes a sum of nearest-neighbor pairs. It is 
known that, at half filling, the ground state of Eq. (1) is in¬ 
sulating with an AF order for any positive Ult. However, it is 
important to consider a non-magnetic Mott transition without 
explicitly introducing magnetic orders, because the essence 
of Mott transitions is independent of magnetism; for instance, 
Mott transitions take place in spinless boson systems of ultra¬ 
cold atomic gases.Actually, nonmagnetic Mott transitions 
were found in the normal and J-wave pairing branches in the 
Hubbard model Eq. (1) at f/ ~ W {W\ band width), and their 
properties were studied using a few methods.Furthermore, 
as we will see later, the AF state also undergoes a crossover 
from a Slater-type (or band) insulator to a Mott insulator at 
U ~ IT. In this work, we mostly treat repulsive cases, but also 
study attractive cases {Ult < 0), in which the ground state is 
5*-wave SC for any Ujt, but a normal state, as an excited state, 
undergoes a spin-gap transition. 

To this model, we apply a variational Monte Carlo (VMC) 
method,^’which yields reliable results for any correlation 
strength UIt and doping rate 6=1- N^IN^ (Nq. number of 
electrons, N^: number of sites) in that the local correlation fac¬ 
tors are exactly treated. As a many-body wave function 'F, the 
Jastrow type is useful: 'F = Here, O indicates a Hartree- 
Fock-type one-body state, and ^ is a product of many-body 
projection factors. In previous VMC studies,Mott 
transitions were found and the aspects of Mott physics were 
revealed using many-body states 'F = or similar wave 

functions. Note that a D-H binding factor Pq plays an essen¬ 
tial role for describing Mott physics.^^ Details will be given 
in the next subsection. 

To accurately compute variational expectation values of 'F, 
we use a variational Monte Carlo method in which the vari¬ 


2.2 Trial wave functions for repulsive cases 

In this subsection, we mention a conventional part of wave 
functions used for systems without a Peierls phase A in pre¬ 
vious studies.For systems with A, we additionally need 
a phase factor Pq, as we will discuss in Sec. 2.3. We study 
Mott transitions and Mott physics in a correlated J-wave sin¬ 
glet pairing state ('Fj = P^d) and a projected Fermi sea 
('Fn = ^^f) as a normal state, and a crossover between a 
band insulator and a Mott insulator for a projected AF wave 
function ('Faf = ^^af)- 

First, we explain the one-body part O. The BCS function 
with dx 2 _y 2 - (uniform s-) wave gap is written as. 







1 |vac). 

(2) 

^<i(5)(k) 


(3) 

- 0^ + |Ad(i)(k)P 


where = -2t(coskx-\-cosky), Ad(k) = A(coskjc-cosk 3 ;) and 
A^ = A. { and A are variational parameters which coincide 
with chemical potential and singlet pairing gap, respectively, 
in the limit of Ult ^ 0. The Fermi sea (FS) is given by Op = 
nkeFS,o- ^ mean-field-type AF state by 


<I>AF = Y\ + sgn(a-)/3k4+Qcr] l^ac), (4) 

keFS,o- 


where Aaf is a variational parameter corresponding to the 
mean-field AF gap, and sgn((T) = +1 according to cr or 

Next, we explain the many-body part P = PqPq. The most 
fundamental onsite (Gutzwilier) projector,^ 

^^0=(6) 


1 

1 - (+)— 

^^k 

2 

^ a/ 

^k ^ ^AF 7 


controls the density of doublons, d = {T(mi)l{UN^)\ as the pa¬ 
rameter g decreases, d decreases. Correspondingly, the range 
of g is 0 < g < 1 (1 < g < oo) for UIt > f) {UIt < U). 
For g = 0, doublons are completely excluded. For Mott phys¬ 
ics, another correlation factor Pq is crucial, which controls 
the binding between a doublon and a holon^’^ and is written 
explicitly as, 

'Pq = n^i-GA (7) 

j 
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Qj = Uddj no hj+j) + l^hhj no -dj^r), (8) 

T T 

where dj = hj = (1 - - nji) and r runs 

nearest-neighbor sites of the site j. To treat doped (asymmet¬ 
ric) cases, we distinguish the contribution of a D-to-H con¬ 
figuration from that of a H-to-D one,^^ which are controlled 
by the parameters jid and jUh, respectively. It was repeatedly 
shown that this type of short-range D-H binding factors cap¬ 
ture the essence of Mott transition and Mott physics,except 
for the Drude weight."^ Although the form of Eq. (8) is slightly 
different from what was used by Millis and Coppersmith,"^ the 
properties of the two are essentially the same. 



Fig. 1. (Color online) Schematic figure of an exchange process in Heisen¬ 
berg model and attached phase factors. 


the AF Heisenberg model (/ > 0), 


"Hneis 


(ij) 





( 12 ) 


2.3 Drude weight and phase factor Vq 

According to Kohn,^ in calculating the Drude weight, we 
add a vector potential A as a Peierls phase to TTkin in Eq.(l): 

-HCA) = -f E + H.c.) + (9) 

<ij>,cr 


Assuming A = Ax, the Drude weight in v direction is given 

by, 


D = 


d^E(A) 


dA^ 


A=0 


( 10 ) 


where E{A) = {AYH{A)\A) with |A) being the normalized 
ground state of El {A). Actually in this study, we obtain D by 
calculating [E{A) - £’(0)]/2A^ for small A’s with A < tt/L (L: 
linear dimension of the system). If the system is insulating 
[metallic], the relation E{A) = Eif)) [E{A) oc A^] holds, and 
D = 0 [D > 0]. 

In constructing a trial state for El {A), we should notice 
that the matrix elements of Eq. (9) with respect to real-space 
configurations are complex. Because, in general, eigenvectors 
of an Hermitian matrix are essentially complex, a trial wave 
function should be complex. To recognize the importance of 
imaginary part in the wave function, we apply Hellmann- 
Feynman’s theorem in the second order to Eq.(lO), 


D = <0| 


cfidiA) 


dA^ 


| 0 > + 2 < 0 | 


dldiA) 


A=0 


dA 


A|A) 
A=0 dA' 


A=0 


<0|rE(4cr^'+- + H.c.)|0) 


-2i{0\tJ]iclci^,a-n.c.) A|A> 


A=0 


( 11 ) 


where |0) is the normalized ground state of El(0). Here we 
assume that the ground state of El(0) does not have a param¬ 
agnetic current, namely, d{A\El(A)\A}/dA\A=o = 0. The first 
term of Eq. (11) is the absolute value of kinetic energy in 
X direction, which is independent of A. On the other hand, 
the second term depends on A and includes imaginary unit i. 
When the system becomes insulating (D = 0) with a finite 
kinetic energy, the first term must be cancelled out by the sec¬ 
ond term. Because D is real, |A) must have a finite imaginary 
part. Thus, it is natural that the ordinary wave functions (^O) 
discussed in Sec. 2.2, including what was used by Millis and 
Coppersmith, exhibit metallic behavior for the Drude weight 
(D > 0) for U/t < oo, because these functions are real. 

Now, let us construct a trial wave function suitable for a 
Mott insulating state with a finite A. We start with considering 


which is the effective model in a Mott insulating regime (U > 
Uc, S = 0) of the Hubbard model. If a finite A is applied— 
in this case the original Hamiltonian becomes Eq. (9)—, the 
second-order virtual hopping processes in the strong-coupling 
expansion illustrated in Fig. 1 does not yield a phase, because 
the two hopping processes occur sequentially, and phase fac¬ 
tors cancel out each other: 

= cJ^cj-^A^cn. (13) 

As a result, Eq. (12) is invariant irrespective of whether A is 
null or finite, and the matrix elements remain real. Thus, a 
phase factor is needless in a wave function for the Heisenberg 
model.^^ 


—E ^ 4" “t" 4" 4" — 

(a) (b) (c) 


Fig. 2. (Color online) Illustration of hopping in ‘H(A) in Mott-insulating 
region. A phase [e~^^] is added when an (t-spin) electron hops to the left 
(a) [right (c)] nearest-neighbor site from the original configuration (b). 


On the other hand, in the corresponding virtual processes 
in the Hubbard model Eq. (9) with a large U/t —precisely, 
the two hopping processes in which a doublon is created and 
annihilated—, the two hoppings do not necessarily occur se¬ 
quentially. Namely, we have to consider the two hoppings 
are mutually independent, although the two events occur in a 
paired manner (if not, the phase becomes metallic). Thus, we 
have to eliminate a phase added in a hopping process shown 
in Fig. 2 within each process. To this end, we introduce into 
the wave function a configuration-dependent phase factor. 


Pe = exp 


W ^ dj (hj+i hj- 1) 

j 


(14) 


where j indicates the coordinate in x direction, and ^ is a var¬ 
iational parameter. If the relation ^ = A is satisfied, the total 
phase factor vanishes in each hopping process. 

This is also viewed from energetics in Eq. (10). Assuming 
the trial wave function is real, kinetic energy increases propor¬ 
tionally to -Re(f^-^^) when an ('[-spin) electron in Fig. 2(b) 
hops to the right or left site. Because an insulating state sat¬ 
isfies E(A) = E(0), this phase has to be cancelled by another 
phase factor such as Eq. (14). 
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If one applies 'F = to 7-f(A) with A = 0, the 

phase parameter 6 in Pg is optimized at ^ = 0, and Pq^q^g 
is reduced to PqPq discussed in Sec. 2.2. Thus, the results 
of quantities other than D obtained in the previous studies are 
not modified by introducing Pq. 

Finally, we discuss some points related to (1) The form 
of Pq in Eq. (14) is the most fundamental one. It is possi¬ 
ble to construct an improved trial state that takes account of 
wider electron configurations. We will return to this point in 
Sec. 3.2. (2) Pq, which is configuration dependent, is concep¬ 
tually different from position-dependent phase factors used in 
different contexts.On the other hand, Millis and Copper¬ 
smith discussed the importance of configuration-dependent 
phase factors for insulating behavior in their original paper,^ 
but they did not provide a tractable form of wave functions. 
(3) In general, this type of configuration-dependent phase fac¬ 
tor seems crucial for constructing current-carrying states in 
regimes of Mott physics in Hubbard-type models which al¬ 
low double occupation. We return to this point in Sec. 6. 

2.4 Drude and SC weights in variation theory 

The expression of the Drude weight D in Eq. (10) is iden¬ 
tical to that of SC (Meissner) weight /)§, but D and are 
different entities, in general. Scalapino, White and Zhang 
(SWZ)^ proposed a criterion of distinguishing them in using 
Eq. (10) by considering the possibility of level crossings at 
small A’s: D is given by the A derivative of the ground-state 
energy at A = 0 (adiabatic derivative), while is given by 
the A ^ 0 limiting value of A derivative of the ground-state 
energy at finite A’s (envelope derivative). 

However, this criterion has subtle points,^^ and is not di¬ 
rectly applicable to variation theory, because, generally, a trial 
function is not assumed to undergo frequent level crossings 
for A ~ 0. If one sets a normal (SC) state as a trial state at 
A = 0, the state remains normal (SC) at a finite A in most 
cases. Therefore, henceforth, we simply regard the quantity 
obtained using Eq. (10) as D (Ds), if a trial state is incoher¬ 
ent (coherent) such as 'Fn and 'Faf ('Fj and 'F^). We consider 
D = for coherent states, and 0^=0 for incoherent states, 
following common relations.^ 

3. d-wa\e pairing and normal states for t/ > 0 

In Sec. 3.1, we discuss the results of J-wave pairing state 
('Fj) and normal state ('Fn) with P = PqPqPq. It is known 
that 'Fj ('Fn) exhibits a first-order Mott transition at Udt ^ 

6.5 (8.6)^^’^^ in quantities such as d, ^(k) and A^(q). We 
mainly display the results of 'Fj, because the behavior of D 
and Mott transitions is qualitatively similar between 'Fj and 
'Fn- In Sec. 3.2, we improve the phase factor upon Pq. 

3.1 Plain phase factor Pq 

First, we show how energy is reduced by introducing the 
phase parameter 6. In Fig. 3, we plot total energy per site E 
for a few values of small A and 6. Here, all variational param¬ 
eters other than 6 are optimized. Shown in the panels (a) and 
(b) are the cases of weak correlations, Ujt = {) and 4, respec¬ 
tively, in which the systems are metallic (fJ < Uf), as men¬ 
tioned. Let the optimized 6 be ^min- We find ^min is situated at 
or in the very vicinity of zero for any values of A and 6. As 
a result, E{A) basically becomes an increasing function of A. 
On the other hand, in a strongly correlated regime {U > Uf) 



0 0.1 0.2 0.3 0.4 0.5 

Ny=A X LUn 

Fig. 4. (Color online) Increment in total energy when a vector potential 
A is applied is shown as a function of A. Cases of four values of Ult are 
compared for d = 0 (half filling) and 6 = 0.1389. Udt ~ 6.5. 



3.0278 

1.0833 

3.1389 

.1944 

3.25 


Fig. 5. (Color online) The optimized phase parameter 6 is plotted as a func¬ 
tion of Up for several doping rates. The arrow indicates the Mott transition 
point at half filling. 


[Figs. 3(c) and (d)], O^m shifts to an appreciably large value as 
A increases [see the arrows in Fig. 3(d)]. However, the value 
of E{A) at ^ = ^min changes only very slightly when A is var¬ 
ied. This situation is summed up in Fig. 4, where we plot the 
energy increment when a small A is applied. Here, all the var¬ 
iational parameters including 6 are optimized. For Ult = 0 
and a large L, E{A) - E(f)) is expanded with respect to A in a 
quadratic function as, 

£o(A)-£o( 0 ) = A 2 r(T\ f dkcosk, + ---. ( 15 ) 

\27r/ JkekF 

First, we consider the half-filled case (solid symbols in 
Fig. 4). We find this quadratic behavior still continues in the 
metallic regime (U/t = 4). On the other hand, in the insu¬ 
lating regime (U/t = 12, 16), E(A) becomes almost constant 
[= £”(0)], as mentioned. Thus, the behavior of E(A) - E(0) 
qualitatively changes through Udt. To pursue it, we show, 
in Fig. 5, the optimized phase parameter ^ as a function of 
U/t. For 6 = Q, 6 exhibits a discontinuity at f/ = t/c ~ 6.5 l 
W e confirmed that this value of Udt precisely coincides with 
that previously estimated from other quantities.Note that, 
for U > Uc, 0 approaches A, and increases as L increases as 
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N{=exLI2n) Af(=6>XZ,/2Tc) 

Fig. 3. (Color online) Total energy per site is plotted as a function of the phase parameter 6 for the J-wave pairing wave function In each panel, data 
for a couple of doping rates and vector potentials are shown; among the four panels, the value of U It is different. We scale 6io N = 6x Ljln = 6x I2l2n in 
abscissa. The arrows in (d) indicate the optimized values of 6 (0min) for d = 0 and the three values of A. 


L=\2 

0.8 


0.6 

Q 

0.4 


0.2 


0 

0.8 

0.6 

Q 
0.4 

0.2 

0 

Fig. 6. (Color online) In (a) and (b), the Drude weight is plotted as a function of Ult for the d-wave pairing and normal states, respectively. The symbols are 
common. In each panel, data for several d’s are plotted. The arrows indicate the Mott transition points at half hlling. For the normal state, Udt ~ 8.7.^^ In (c) 
and (d), the same quantity is shown as a function of doping rate for the same states, respectively. In each panel, finite-C//r VMC data (L = 12) are shown with 
symbols; Solid symbols and solid lines denote the results of and open symbols and dashed lines those of discussed in Sec. 3.2. For comparison, the 
behavior of D for [//r = 0 (L = oo) is plotted with a black solid line. 
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shown in the left inset of Fig. 7, in contrast to the behavior 
for U < Uc, where 6 decreases with L. It demonstrates that 
the mechanism assumed when the phase factor Pq was intro¬ 
duced [Eq. (14)] certainly works for U >Uc. 

We turn to doped cases. As shown in Fig. 4, E{A) - £’(0) 
continues to be a quadratic function of A even for sufficiently 
large values of UIt. In Fig. 5, no discontinuity in the behav¬ 
ior of optimized 6 is found even for the smallest doped case 
{6 = 0.0278). Thus, the Mott transition vanishes on doping, as 
expected. As 6 increases, the optimized value of 6 decreases 
from A. This is mainly because an isolated (untied to dou- 
blon) doped hole adds a phase A or -A during the hopping, 
but Pe does not compensate for it. In general, it is probably 
impossible to make a phase factor which compensates for the 
phase generated by a free carrier. We may say this is the ori¬ 
gin of conduction. Another reason for ^ < A is interplay of 
a doublon and multiple holons; we will take up this effect in 
Sec. 3.2. 



Ult 


d=0 


L 

—*—10 

-^12 

-▼-14 

16 


Fig. 7. (Color online) Magnification of Drude weight for d = 0 in J-wave 
pairing state near Mott transition points Udt ~ 6.5 for several values of L. In 
the right inset, the system-size dependence of D is shown for four values of 
Ujt in the Mott-insulating regime. We estimate D with A = InjLxO.l. In the 
left inset, the system-size dependence of optimized 6 is plotted versus UIt. 


Using E{A) obtained above, we estimate the Drude weight 
through Eq. (10). In Fig. 6, we plot D for the J-wave pairing 
and normal states (solid symbols). To begin with, we consider 
Ult dependence [panels (a) and (b)] at half filling. For Ult = 
0, as derived from Eq. (15), the Drude weight becomes, 

, 16 , 

where is the kinetic energy in x direction. Thus, we have 
D = 8/71^ at d = 0 for the square lattice. As Ult increases, 
D rapidly decreases and almost vanishes dXU = Uc in both 
states as shown by arrows in Figs. 6(a) and 6(b). In Fig. 7, 
the behavior of D near the Mott transition point is compared 
among different L’s. It reveals D exhibits a discontinuity at 
U = Uc for large L’s (L > 12).^"^ Note that D decreases as 
L increases for U > U^ in contrast to the feature for U < 
Uc. In the right inset of Fig. 7, D is plotted as a function of 
IIL for some values of Ult for U > Uc. Although accurate 
extrapolation is not easy owing to the statistical fluctuation for 
large L's, D probably vanishes for L ^ oo. The value of Udt 
and the features of a first-order transition precisely coincide 
with those obtained in other quantities. Thus, we could show 


for and 'Fn that a conductor and a Mott insulator can be 
clearly distinguished by the Drude weight, using a variation 
theory. 



Fig. 8. (Color online) Contribution of the second term in Eq. (11), namely, 
the difference between the absolute value of kinetic energy in x direction 
and D is plotted for as a function of Ult for several doping rates. For 
comparison, \E^^\lt for d = 0 is shown with a dash-dotted line. The arrow 
indicates the Mott transition point at half filling. 


Now, we look at the effectiveness of Pq with respect to Ult 
and d. As discussed, without the imaginary part in 'F, D is 
reduced to |L^^|/f [the first term in Eq. (11)], which is pro¬ 
portional iotlU fovU It ^ oo. The contribution of the second 
term in Eq. (11), which is the direct effect of Pq, is given 
by |L^.^|/f - D. In Fig. Ult dependence of this quantity is 
shown for some d’s. At half filling, this quantity is very small 
compared with \E^JIt for U < U^ whereas the two quantities 
approach each other for U > Uc. As S increases, the contri¬ 
bution of \E^JIt - D for U > Uc rapidly becomes weak, in 
contrast to the increase in |L^.^|/L Consequently, D increases 
as S increases, as we see next. 

Shown in Figs. 6(c) and 6(d) is the doping dependence of 
D. For Ult = 0 (black solid line), the Drude weight starts 
from 8 171^ at d = 0 and monotonically and slowly decreases 
as S increases, because Eq. (16) holds. Eor U < Uc, this weak 
dependence on S continues as seen for Ult = 4. Here, the 
alternate behavior for a small Ult stems from the boundary 
conditions we use (a finite-size effect), and is irrelevant. The 
behavior of D abruptly changes around U = Uc. Eor U > 
Uc, D increases linearly as S increases. Difference is slight 
between the normal and J-wave SC states. 

Linear behavior of D(6) was observed for the t-J model,^^ 
and therefore this behavior is characteristic of doped Mott in¬ 
sulators. Also in strongly correlated Hubbard models, the be¬ 
havior of D(d), namely, D oc for d ^ 0 has been studied 
for long, mainly on the basis of exact diagonalization.^^“^^ 
The main concern of these studies was whether the exponent 
a is 1^^’^^ or 2?^ On the other hand, it is well-known that 
linear behavior of D{6) was found in the cuprate supercon¬ 
ductors first by yuSR experiments (so-called Uemura plot),^^ 
which gave strong evidence that the cuprate SC’s are doped 
Mott insulators. A recent experiment"^^ showed SC weight ps 
looks somewhat convex, which is similar to for U > Uc in 
Eig. 6(c). Although the relationship of the present results to 
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experiments has to be deliberately analyzed, we believe both 
capture a typical feature of Mott physics. 

3.2 Improved phase factor 
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Fig. 9. (Color online) Difference of attached phases is illustrated between 
hopping processes from (a) to (b) and from (c) to (b), which make a doublon 
with two neighboring holons in x and y directions [(b)]. 


Average distance between a doublon and a holon becomes 
appreciably larger than the nearest-neighbor-site distance for 
U ^ However, when we construct we only take ac¬ 
count of the contribution of nearest-neighbor D-H pairs. In 
addition, when holes are doped, the interplay of a doublon 
and multiple holons becomes effective. As shown in Fig. 9(b), 
a doublon is frequently next to two holons in v and y direc¬ 
tions simultaneously for d > 0. This configuration is real¬ 
ized by a single hopping from a configuration in Fig. 9(a) or 
Fig. 9(c). Although these two hoppings occur with the same 
possibility if A is sufficiently small, added phases are different 
between the two, namely, e~^^ (e^) is added in the former (lat¬ 
ter) hopping. On the other hand, the phase factor Pg attaches 
the same counter-phase 6 to the two hopping processes. Thus, 
even without considering the effect of isolated holons, ^min de¬ 
viates from A SLsUlt becomes small or S increases. It is better 
to treat a D-H pair in x direction independently of a pair in 
y-direction. 
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Fig. 10. (Color online) We illustrate the way of assigning phase factors 
(j = 1,2, • • • , 8) to local electron configurations in the improved phase factor, 
p 0 . Details are explained in the text. 


On the basis of this argument, we substitute an improved 
projection factor Pq for Pq. In Pq, we allow for multiple- 
holon effects as well as extended D-H binding effects up to 
the two-step distance with eight phase parameters (Oi-O^) to 


Table I. Improvement of total energy for the normal state (L = 12) by 
Po on pQ is summarized, when a Peierls phase A [= 0.5L/(27r)] is applied for 
some values of Ujt and 6 . The second, third, and fourth rows indicate the total 
energies of the given wave functions. The values in brackets indicate errors in 
the last digits. The last low shows the ratio [E(^)-E{Pe^)\l[E(^)-E{Pe^)\. 


Up, 5 

4, 0.0 

12, 0.0 

12, 0.1389 

16, 0.0 


-0.77484(7) 

-0.2538(2) 

-0.5779(1) 

-0.2029(2) 


-0.77551(7) 

-0.2608(9) 

-0.5796(1) 

-0.2092(1) 


-0.77599(7) 

-0.2610(2) 

-0.5808(1) 

-0.2091(2) 

Ratio 

0.58 

0.97 

0.59 

1.00 


be optimized. Because the operator representation of Pq is 
complicated, we explain the function of Pq with illustrations 
in Fig. 10. Formally, Pq is represented as 

£ 

where ^ is a site index. For example, if a doublon at site i 
(orange circle in Fig. 10) has a nearest-neighbor holon as 
in Fig. 10(a), e^^^ is assigned. If a doublon has no nearest- 
neighbor holon but a second-neighbor holon as in Fig. 10(d), 
e^^^ is assigned. If a doublon has no first- and second-neighbor 
holon, as well as if site i is empty or singly occupied, we as¬ 
sign e^^. In Fig. 10, we only explain D-to-H factors (ffi-^ 4 ), 
but we also consider H-to-D factors (Os-O^) in the same way. 
In Fig. 10, we only show configurations in which D is located 
on the left to H; when the positions of D and H are exchanged, 
an assigned phase factor should be e~^^J. 

The Drude weight calculated with Pg^ is shown in Fig. 
6 with open symbols. As compared to the results of simple 
(solid symbols), the magnitude of D becomes small for 
any values of Ult and 6 . For a quantitative analysis, we com¬ 
pare, in Table I, the total energy for a small A among three 
wave functions, namely, 'Tn without a phase factor, Pf^^ 
and Pe^^ for typical values of Ult and 6 . In the last row, 
we show the ratio of improvement by Pf^I to the improve¬ 
ment by Pe^. Although the statistical fluctuation is large, we 
can grasp a tendency. The improvement by Pq^ is notable for 
an intermediate U/t or sl finite S; not a small portion [~ 40% 
for (U/t, S)=(4, 0.0) and (12, 0.1389)] of the improvement is 
achieved by the newly introduced part in Pe^. On the other 
hand, for d = 0 and U > Uc, the difference between E(P 0 ^) 
and E(p 0 ^^) is negligible. This result is what we expected. 

In summary, the essence of Mott transitions as to the Drude 
weight is captured by the simple phase factor Pq, but quan¬ 
titative improvement is possible for intermediate correlation 
strengths or finite doping rates by introducing refined phase 
factors such as Pg. 

4. Antiferromagnetic state 

In this section, we study the Drude weight of AF states. In 
contrast to the normal and J-wave pairing states, which are 
insulating only for U > Uc Sit half filling, the AF states are 
insulating for any positive value of Ujt for d = 0, because 
the nesting condition is completely satisfied for hypercubic 
lattices (cf. Fig. 13). Thus, the nature of AF states as insulators 
changes from a band (Slater) insulator to a Mott insulator as 
Ult increases. As we will see shortly, AF states become the 
touchstone of applicability of the phase factor Pq. 

Let us start with the results of 'Tap obtained by VMC. In 
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Fig. 11. (Color online) The Drude weight of T^af at half filling is shown as 
a function of UIt. Solid (open) circles indicate VMC results of using Pq (Pq) 
as a phase factor. For comparison, \E^^J/t is shown with open triangles. 



Fig. 12. (Color online) Drude weight of T^af calculated using VMC as a 
function of doping rate. Solid (open) symbols indicate the results of using Pg 
(Pg) as a phase factor. For large d’s (e.g. d > 0.15 for t//r = 12), the state 
becomes normal. 


Fig. 11, we plot Ujt dependence of D at half filling. Similarly 
to 'Fn and 'Fj IFigs. 6(a) and 6(b)], D starts from D = S/tt^ for 
f//f = 0 and decreases, as Ult increases, rapidly but smoothly 
owing to absence of a transition. First, we discuss the strongly 
correlated regime. For UIt > 10, D becomes very small, com¬ 
pared with indicating that Pq appropriately works also 

for 'Faf in this regime. Doping-rate dependence of the Drude 
weight shown in Fig. 12 is almost linear, and similar to those 
for 'Fn and 'Fj IFigs. 6(c) and 6(d)]. Next, we consider the 
weakly correlated regime {U < W), where 'Faf is a band in¬ 
sulator. Being an insulator is confirmed by the behavior of the 
momentum distribution function and charge density structure 
factor shown in Fig. 13. Namely, there is no Fermi surface 
Ino discontinuity in ^(k)], and a charge gap opens tA^(q) cx 
with Of > 1 for |q| ^ 0]. Therefore, D should vanish even for 
U < W. However, the Drude weight calculated with 'Faf ex¬ 
hibits large finite values for U <W. This indicates that 'Faf is 
lacking in some important element for describing D in a band 
insulator. 

It is useful to consider this point with a mean-field theory. 



Fig. 13. (Color online) (a) Momentum distribution function and (b) charge 
density structure factor of'P af (solid symbols) and 'Fn (open symbols) calcu¬ 
lated using VMC are plotted along the path (0,0)-(7r, 0)-(7r, 7r)-(0,0) for three 
values of Ult. The Mott transition point for 'Fn is Udt ~ 8.7. 


assuming Ujt is small. Note that the one-body part Oaf is an 
insulating wave function for t//f > 0 at half filling, even with¬ 
out correlation factor !P. Therefore, D has to vanish within the 
mean-field theory, meaning the one-body part has to be com¬ 
plex. Because Oaf is essentially real, we need to modify it in 
this context. Let us look at this point in the light of a pertur¬ 
bation theory. We expand 7T(A) in Eq. (9) with respect to A 
as, 

7T(A)^<K(0)-f A/;,, (18) 


with paramagnetic current 

Jx = -it L(4crCr+xo- - H.C.). (19) 

r,cr 


Assuming A is small in Eq. (18), the one-body wave function 
within the first-order perturbation is given by 


= |0) + a ^ 


{fn\Jx\0) 

Eq — Efn 


|m>. 


( 20 ) 


where |0) is the ground state of the mean-field Hamiltonian 
('Hmf) (equivalent to Oaf), N) represents the excited states 
of and Eq and E^ are the corresponding energies. To 

consider the perturbed term in Eq. (20), we apply a Eourier 
transformation to J/. 

Jx — ^ sin ^x ^^j ^qcr^qo'’ (^1) 

qx qy,o- 
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We may replace Cq^Cqo- with the operators of Bogoliubov 
quasiparticles for Fermi sea, J-wave BCS or AF states, ac¬ 
cording to the choice of mean fields. Because Eq. (21) is al¬ 
ready a diagonal form for the Fermi sea and BCS’s Bogoli¬ 
ubov quasiparticles, {m\Jx\^) vanishes. Consequently, the nor¬ 
mal and J-wave pairing states do not have a first-order cor¬ 
rection in Eq. (20). On the other hand, if we substitute AE- 
Bogoliubov quasiparticles in Eq. (21), we have 

Jx = -2r E sin^x[(Q'q-/3q)(aqcraq(r-aq+Qo-aq+Qcr) 

qgMBZ 

—2(Tafq/5q(^Zq^^2q+Qcr + ^q+Qcr^qcr)j , (22) 

where MBZ indicates the reduced magnetic (AE) Brillouin 
zone. Owing to the interband term ^q+Q^^qo-, {MJx\^) hoes 
not vanish in this case. In summary, for band-insulating states 
such as the AE state, we need to add an imaginary corrective 
term in Eq. (20) to the ordinary one-body part. 



Fig. 14. (Color online) Doping-rate dependence of Drude weight of AF 
state estimated using a mean-field theory for a couple of small values of t//r. 

Within the present VMC formalism, it is not easy to calcu¬ 
late D directly using a wave function like Eq. (20). We need 
to develop a technique to treat multiple determinants. Instead, 
we can calculate D of the AE state within a mean-field the¬ 
ory.^ In Eig. 14, we show doping-rate dependence of D thus 
estimated for small values of Ult. At half filling, vanishing of 
D is realized. As d increases, D should linearly increases for 
small d’s, assuming that the AE state is intrinsically stable. 
However, it is known that doped AE states are unstable to¬ 
ward phase separation for the simple square lattice,^^ so that, 
actually, AE states are not realized for d > 0. Anyway, this 
analysis of AE states reveals that the mechanism of reducing 
D in band insulators like AE states is distinct from that in Mott 
insulators. 

5. Attractive Hubbard model 

In this section, we study the behavior of SC weight and 
Drude weight D in the attractive Hubbard model {Ult < 0), 
when the configuration-dependent phase factor Pq and a re¬ 
fined version is applied. 

Using a canonical transformation on a bipartite lattice, 
one can map the physics of repulsive Hubbard model at half 
filling in magnetic fields to the physics of attractive Hubbard 


Table II. Comparison of total energy [L = 12, A = 0.5 x (27r/L)] of SC 
state among cases with and without Pq and of correlator product state (CPS) 
for attractive interaction. The values in brackets indicate errors in the last 
digits. 



U/t 


6 

-5 

-10 

-20 

no Pe 

0.0278 

-3.13694(5) 

-5.25918(9) 

-9.93121(6) 

0.0833 

-2.99487(6) 

-4.97905(6) 

-9.37431(5) 

0.1389 

-2.84941(5) 

-4.69652(6) 

-8.81596(5) 

0.1944 

-2.70042(5) 

-4.41163(6) 

-8.25631(5) 

n 

0.0278 

-3.13705(6) 

-5.25923(6) 

-9.93126(8) 

0.0833 

-2.99501(5) 

-4.97908(5) 

-9.37431(6) 

0.1389 

-2.84955(5) 

-4.69659(7) 

-8.81601(6) 

0.1944 

-2.70055(4) 

-4.41171(6) 

-8.25635(5) 

CPS 

0.0278 

-3.13717(5) 

-5.25942(5) 

-9.93149(7) 

0.0833 

-2.99507(5) 

-4.97926(7) 

-9.37456(5) 

0.1389 

-2.84965(5) 

-4.69675(7) 

-8.81623(5) 

0.1944 

-2.70065(4) 

-4.41185(6) 

-8.25658(6) 


model.^^ By applying this transformation, the form of Gutz- 
willer projection Pq remains intact (1 < g < oo), but the D-H 
binding projection for f//f > 0 is transformed to 

pQ = (23) 

j 

T T 

for Ujt < 0.^^ Here, = njcr{\ - ^y-cr) and yu (0 < // < 1) is a 
variational parameter which controls the binding between up- 
and down-spin electrons. In a previous paper,^^ it was shown 
using VMC calculations that the normal state 'Fn = 
undergoes a transition between a metallic and a spin-gapped 
phases at lUd/f ~ 9 irrespective of electron density. On the 
other hand, a homogeneous 5--wave singlet pairing state = 
pQ^G^s exhibits a crossover in the SC properties from a BCS 
type to a Bose-Einstein condensation type at \ Uco\lt ~ 8.7. 

The SC weight was calculated using and 'F^ = 

pQ^G^s^^ as functions of Ujt (< 0). Because these wave 
functions are real, resultant Ds’s are substantially a half of the 
absolute values of kinetic energies (see, for instance, Eig. 18 
in Ref.^^). Nevertheless, the results are broadly consistent 
with those of quantum Monte Carlo^^ and DMET^^“^^ calcu¬ 
lations. Thus, the improvement by phase factors is expected 
to be small in contrast to the repulsive case. Here, we check 
the effectiveness of the phase factors for Ujt <{t. 

As a first choice, we consider the same form Pq in 
Eq. (14),^^ and calculate Ds {D) with For a 

large |t/|/f, most electrons form doublons, but Pq does not 
act on the hopping processes concerning isolated doublons. 
Eurthermore, Tq does not cancel the phase +2A of a doublon 
hopping Thus, the effect of Ve is expected to be 

limitative. In Table II, we compare £’/r between the cases with 
and without Pq for some values of 6 and UIt. We find that 
E (equivalently D{) decreases by applying but the decre¬ 
ment is extremely small for any Ult and d (orders of 10“^F 
lO'^O. ns expected. To corroborate the fact that phase factors 
are ineffective in this case, we introduce a refined phase fac¬ 
tor in terms of the correlator product state (CPS).^^’^^ The 
phase factor used here depends on the electron configurations 
of local five-site clusters, and has 4^ phase parameters to be 
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optimized. As shown in lower rows of Table II, even the re¬ 
sults of CPS improve the energy only slightly. In conclusion, 
a configuration-dependent phase factor is not an important el¬ 
ement of Ds for < 0. 

6. Summary and Discussions 

In this study, we considered how to appropriately calcu¬ 
late the Drude and SC weights in the variation theory for the 
Hubbard model. We argued that existence of the imaginary 
part is indispensable in the wave function for suppressing D 
(and Ds). In strongly correlated regimes where Mott physics 
prevails, a phase correlation factor YPq in Eq. (14)], which de¬ 
pends on the local electron configuration, works successfully 
to estimate D and in normal and SC states, respectively. 
Thereby, at half filling, D vanishes for U > Uc, where Udt 
is the Mott transition points previously determined by other 
quantities such as doublon density and charge susceptibility 
[Figs. 6(a) and 6(b)]. Namely, Mott transitions are correctly 
described using D. Thus, the long-standing problem of Millis 
and Coppersmith^ was resolved. 

Doping-rate dependence of D and for f/ > f/c is linear 
for d ^ 0 and widely an increasing function of d with some 
convexity [Figs. 6(c) and 6(d)]. This behavior is consistent 
with what is observed experimentally for cuprate SC’s.^^’^^ 
On the other hand for weak correlations {U < Uc), D is finite 
at half filling, and only weakly dependent on S especially for 
A. 

As for AF states, which is insulating for any positive Ult 
at half filling in the hypercubic lattices, Vq is still effective 
for U > W, and D almost vanishes. However, for f/ < W 
[band-insulating (Slater) regime], D estimated using 
becomes finite and smoothly approaching the value of nonin¬ 
teracting (U = 0) limit di^ Ult decreases. To correct 

this flaw, a perturbation theory with respect to A is useful. We 
revealed that the one-body part of'T should have a finite imag¬ 
inary part, which is introduced by the first-order perturbation 
[Eq. (20)]. This contribution survives for the AF states, but 
vanishes for the normal and SC states. Thus, it can be shown 
D vanishes for Oaf at half filling within the mean-field the¬ 
ory.^ In a band insulator like 'Taf with U Uc, D vanishes 
through a mechanism different from the doublon-holon bind¬ 
ing effect in a Mott insulator. 

For the attractive Hubbard model, we checked the effect of 
configuration-dependent phase factors in calculating D (D^) 
in normal (^--wave SC) states. The effect of phase factors ex¬ 
ists, but is considerably small as compared to the repulsive 
model of a large UIt. 

In the remainder, we make a few discussions related to the 
present subject. 

(i) The Drude (or SC) weight is an important quantity, but, 
technically, D does not seem a very suitable measure to dis¬ 
tinguish a conductor from an insulator at least in the varia¬ 
tion theory, because D requires fine tuning of the imaginary 
part of the trial states according to the situations of individual 
states. As an alternative measure, localization length T, which 
estimates the degree of insulating and diverges in conductive 
phases,^^ is more convenient in the sense that special tuning 
is not necessary for the trial states in appropriately computing 

(ii) As mentioned, it was shown by introducing a 
configuration-dependent phase factor like Pq that a staggered 


flux (or J-density wave) state, in which a local circular cur¬ 
rent flows in each plaquette of the square lattice alternately, is 
considerably stabilized as compared to the corresponding nor¬ 
mal state in a strongly correlated Hubbard model. Without 
the phase factor, we have no energy reduction. The physics 
in this case is similar to the present one; the vector potential 
of the staggered flux is almost (for d = 0) or partially (for 
d > 0) cancelled by Pq. Similar results were also obtained 
for 2 i d-p model^^ and a Bose Hubbard model^^ with strong 
correlations. It is probable that a phase factor like Pq is gen¬ 
erally necessary for describing current-carrying states in the 
regime where Mott physics is relevant. Inversely, if we can¬ 
not find an effective phase factor for a current-carrying state, 
it is probably not stabilized. 

(iii) In this paper, we chiefly treated a fundamental phase 
factor Pe, which captures the essence of Mott physics. For 
quantitative improvement, it is intriguing to adopt refined 
techniques such as CPS discussed in Sec. 5. 
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